#install.packages("plotly")
#install.packages("ggplot2")
library(plotly)
## Loading required package: ggplot2
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
library(ggplot2)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──
## ✔ tibble 3.2.1 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ✔ purrr 0.3.5
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks plotly::filter(), stats::filter()
## ✖ dplyr::lag() masks stats::lag()
1. Create surface and cutting plane
compose1 <- function(f, c) {
g <- function(y) f(c, y)
return(g)
}
compose2 <- function(f, c) {
g <- function(x) f(x, c)
return(g)
}
composefh <- function(f, h) {
g <- function(x) f(x, h(x))
return(g)
}
get_curve <- function(M, v, f) {
x0 <- M[1]
y0 <- M[2]
a <- v[1]
b <- v[2]
if (a == 0 & b != 0) {
g <- compose1(f, x0)
id <- 1
} else if (a != 0 & b == 0) {
g <- compose2(f, y0)
id <- 2
} else {
h <- function(x) y0 + b * (x - x0) / a
g <- composefh(f, h(x))
id <- 3
}
return(list(g = g, id = id))
}
get_plane <- function(M, v, id, xx, yy, zz) {
x0 <- M[1]
y0 <- M[2]
a <- v[1]
b <- v[2]
if (id == 1) {
YZ <- expand.grid(yy, zz)
X <- matrix(rep(x0, nrow(YZ)), nrow = nrow(YZ), ncol = ncol(YZ))
Y <- matrix(YZ[,1], nrow = nrow(YZ), ncol = ncol(YZ))
Z <- matrix(YZ[,2], nrow = nrow(YZ), ncol = ncol(YZ))
} else if (id == 2) {
XZ <- expand.grid(xx, zz)
X <- matrix(XZ[,1], nrow = nrow(XZ), ncol = ncol(XZ))
Y <- matrix(rep(y0, nrow(XZ)), nrow = nrow(XZ), ncol = ncol(XZ))
Z <- matrix(XZ[,2], nrow = nrow(XZ), ncol = ncol(XZ))
} else if (id == 3) {
XZ <- expand.grid(xx, zz)
X <- matrix(XZ[,1], nrow = nrow(XZ), ncol = ncol(XZ))
Y <- matrix(y0 + b * (X - x0) / a, nrow = nrow(XZ), ncol = ncol(XZ))
Z <- matrix(XZ[,2], nrow = nrow(XZ), ncol = ncol(XZ))
} else {
X <- NULL
Y <- NULL
Z <- NULL
}
return(list(X = X, Y = Y, Z = Z))
}
f_surf <- function(x, y) {
2 + 0.1 * x + 0.3 * y + 0.1 * x^2 - 0.3 * x * y + 0.2 * y^2
}
xx <- seq(-3, 3, length.out = 300)
yy <- seq(-3, 3, length.out = 300)
x <- outer(xx, yy, FUN = function(a, b) a)
y <- outer(xx, yy, FUN = function(a, b) b)
z <- outer(xx, yy, FUN = f_surf)
zz <- seq(1, 5, length.out = 100)
M <- c(1, 1, 0)
v <- c(1, 1, 0)
curve_data <- get_curve(M, v, f_surf)
g <- curve_data$g
id <- curve_data$id
plane_data <- get_plane(M, v, id, xx, yy, zz)
X <- plane_data$X
Y <- plane_data$Y
Z <- plane_data$Z
2. Plot
figure <- plot_ly() %>%
add_surface(x = x, y = y, z = z, colorscale = "Viridis", showscale = FALSE) %>%
add_surface(x = xx, y = yy, z = Z, colorscale = list(c(0, 1), c("rgb(254, 254, 254)", "rgb(254, 254, 254)")), showscale = FALSE, opacity = 0.65) %>%
add_trace(x = X[, 1], y = Y[, 1], z = Z[, 1], type = "scatter3d", mode = "lines", line = list(color = "rgba(50, 50, 50, 0.5)", width = 2)) %>%
layout(scene = list(xaxis = list(nticks = 5, range = c(-2, 2)),
yaxis = list(nticks = 5, range = c(-2, 2)),
zaxis = list(nticks = 5, range = c(1, 5), tickmode = "linear", tick0 = 1, dtick = 1), showbackground = FALSE),
margin = list(r = 20, l = 10, b = 10, t = 10),
font = list(family = "Times New Roman"))
figure